Low temperature solution of the Sherrington-Kirkpatrick model 
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We propose a simple scaling ansatz for the full replica symmetry breaking solution of the 
Sherrington-Kirkpatrick model in the low energy sector. This solution is shown to become ex- 
act in the limit a; ^ 0, fix — > oo of the Parisi replica symmetry breaking scheme parameter x. The 
distribution function P{x,y) of the frozen fields y has been known to develop a linear gap at zero 
temperature. We integrate the scaling equations to find an exact numerical value for the slope of 
the gap dP{x,y)/dy\y^o — 0.3014046... We also use the scaling solution to devise an inexpensive 
numerical procedure for computing finite timescale {x = 1) quantities. The entropy, the zero field 
cooled susceptibility and the local field distribution function are computed in the low temperature 
limit with high precision, barely achievable by currently available methods. 
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The field of spin glass pliysics has been actively stud- 
ied for tlie last thirty years. Much attention has been 
devoted to the Sherrington-Kirkpatrick (SK) model 0|- 
a mean field analog of realistic spin glasses. The mean 
field treatment of this model, although exact, is highly 
nontrivial 2j. In its static formulation, the mean field 
description of the glassy phase involves infinitely many 
steps of replica symmetry breaking (RSB), compactly 
written as a set of integro-differential equations 3] for 
the spin glass order parameter q(x) and some auxiliary 
functions jjj m{x,y) and P{x,y). A dynamical formula- 
tion Q leads to similar results and interprets the variable 
X € [0, 1] as a parametrization of (diverging) timescales of 
the model, with smaller values of x corresponding to the 
longer times a|- The functions m{x,y) and P{x,y) are 
interpreted |7|, respectively, as a local magnetization in 
the presence of a frozen field y and a distribution function 
of the frozen fields, measured at a timescale x. The full 
RSB equations cannot be solved analytically, and even 
solving them numerically is not a trivial task. There were 
attempts to conjecture certain scaling laws 0, which, if 
correct, would grant an exact solution of the problem. 
Those conjectures, however, turned out to be merely a 
good approximation 9]. 

The idea of scaling, nevertheless, deserves attention. 
The spin glass phase is characterized by marginal criti- 
cality, and critical systems generically exhibit some kind 
of universal (scaling) behavior 10^. To be more spe- 
cific we use insights gained from both the static and dy- 
namic viewpoints of the problem. The space of spin glass 
states is ultrametric and can be visualized as a tree of 
states im with the leaves representing quasi-equilibrium 
(pure) states. In the dynamical description the system 
can explore the tree, reaching the states with smaller and 
smaller overlap q, connected to the initial state through 
larger and larger branches. The self-similarity of the tree, 
away from its leaves {x ~ 1) and its root {x <T), is an 
important ingredient which could be related to the scal- 
ing behavior of the functions m(x,y) and P{x,y) in the 
range T <ti x <ti 1. Another condition for universality 
to occur is that the typical frozen fields should play no 



role in the problem. In other words, the scale of disorder 
should be irrelevant. This can only happen in the zero 
temperature limit, when any finite frozen field completely 
polarizes a spin, thus excluding self from contributing to 
the dynamical evolution of the system. Therefore, it is 
the limit of a; — > and Px — > cxd, where one should ex- 
pect to find a scaling solution, which would be valid for 
y 1. In this paper we present such a solution and prove 
that it becomes exact in the above limit, by investigating 
the correction to scaling. We then use the scaling solu- 
tion to solve the Parisi equations on the scale y ^ T and 
X > 0, for the first time directly in the zero temperature 
limit. Various physical quantities, measurable on short 
time scales, are also computed. 

Model. We consider the Sherrington-Kirkpatrick 
model-perhaps the most studied model of spin glass 
physics: 
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where the spins are classical variables = ±1 and the 
bonds are randomly Gaussian quenched with a variance 
J, which is set to unity through the rest of the paper. 
Upon cooling, the ergodicity breaks down at T = 1 and 
the systems freezes in one of the multitude of quasiequi- 
librium states. This phase is referred to as the glassy 
phase and can be described with use of the Parisi RSB 
formalism. A spin on a site i experiences a local field 
J2j Jij ^3 ■ The local fields can be measured with or with- 
out the site i present in the lattice. They are called, ac- 
cordingly, the instantaneous or frozen fields. For T < 1 a 
pseudogap in the distribution of fields opens, becoming 
(asymptotically) linear in the limit ?/ — > 0, — > oo. 

Scaling ansatz. The Parisi infinite RSB scheme equa- 
tions, written in a usual form read: 

q(x) 



m(x, y) 
P{x,y) : 



q{x) 



[m"{x, y) + 2j3xm,{x, y)m'{x, y)] , (2) 
{P"{x, y) - 2(3x [m{x, y)P{x, y)]'} , (3) 



where the dot and prime are derivatives with respect 
to X and y variables, correspondingly. The differential 
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equations are supplemented with the initial conditions: 
m{l,y) = tanh/Sy and P{0,y) = S{y). In principle, these 
equations can be solved iteratively. One can compute m 
and P for a given order parameter q, which, in its turn, 
is computed from m and P: 

q{x) = J dyP{x,y)m^{x,y). (4) 

We introduce new notations: m{x,y) — rh{x,z), 
P{x,y) = {(3x)^^p{x, z), q{x) — 2f]{f3x)^^c{x), where 
z = f3xy. It will become clear later that in the scaling 
regime the functions m, p and c lose their dependence on 
the variable x. We recast Eqs. (|2I3|I using new definitions: 

xm = — c (to" + 2mm') — zrh' , (5) 
xp — c [p" — 2(pto)'] — zp' + p, (6) 

where the dot and prime are now derivatives with re- 
spect to x and z variables, correspondingly. The func- 
tions' arguments arc omitted for compactness. Our scal- 
ing ansatz states that to and p are functions of the 
scaling variable z only, and c is a constant, given by 
c = / dzp{z)(\ — fr?[z)\ as follows from Eq. Q). Here- 
after we refer to Eqs. I|5I6|I . with the left hand sides set 
to zero and c being a constant, as the scaling equations. 

As expected, the ansatz does not respect the boundary 
conditions (where the tree of states is not self-similar), 
but we will demonstrate that it becomes asymptotically 
exact in the scaling regime described above. This is 
precisely where the slopes of the linear gap are formed. 
Thus the scaling ansatz allows us to compute the slope of 
P(l,?/) in the low temperature limit. Because the scal- 
ing equations can be easily integrated (numerically), one 
can obtain the value of the slope with arbitrary preci- 
sion. The initial conditions for the scaling equations fol- 
low from the definition of to and P and from the linearity 
of the gap at large z: 

z = 0, TO = 0, p' = 0; (7) 
|2;|>1, TO = sign(z), p = 7(|z| + 2c). (8) 

With a fixed constant c the function p enters the equa- 
tions linearly, therefore one can set 7=1 when solving for 
c, and then compute the slope as 7 = c[/ dzpr^=\(z)(\ — 
TO^(z))]^^. Up to ten digits of precision we found c = 
0.4108020997, 7 = 0.3010464715. The functions to(z) 
and p(z) are shown in Fig. l(T|l. 

Stability of the scaling solution. To substantiate our 
findings we have to prove the existence of the scaling 
regime. Equations H5I6|I are free of singularities, and be- 
cause the scaling solution differs from the initial condi- 
tions at a; = 1, our scaling [coinciding with Parisi-Toulose 
scaling Q for q{x)] cannot hold for all x. The important 
question is if there is at all such an x where the pro- 
posed scaling solution becomes exact. One can answer 
this question by investigating the stability of the scaling 




FIG. 1: Scaling functions rh{z) and p{z), obtained by integrat- 
ing the Parisi equations [Eqs 15l6ll l under the assumptions of 
the scaling ansatz. The field distribution function is asymp- 
totically linear, p{z) = jIz] + 27c, \z\ S> 1. 

solution, expanding around it to linear order. We look for 
a correction to scaling in the form: drh{x, z) = x^5'm{z), 
Sp{x, z) = x^6m{z) and 6c{x) = x^Sc. This procedure 
leads to a set of inhomogeneous linear differential equa- 
tions 

XSfh = —c [Sin" + 2((5toto)'] — zSfh' + zrh' /c, (9) 
(A - l)Sp = c [5p" - 2{5pm)' - 2{p5m)'] 

— z5p + {zp — p)/c, (10) 

supplemented by two constraints: 

1 = ^1 - -^^ j dz [5p{l - TO^) _ 2pto(5to] , (11) 

0^ J dz [Sp{mf + 2pm 6m] , (12) 

where Sm = Sfii/dc^ dp — dp/dc, and all functions depend 
on the scaling variable z only. The least stable solution 
corresponds to the smallest A. Integrating (numerically) 
the above equations we found A « 5.41, a fairly large 
value, indicating that not only the scaling solution is lo- 
cally stable in the limit a; ^ 0, but also that it serves 
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FIG. 2: Normalized correction to scaling functions Sfh{z) — 
5rh{x, z)/5c{x) and Sp{z) = Sp{x, z)/5c{x), obtained from the 
linear stability analysis Eqs (191121 . Correction to scaling 
terms exhibit quick decay with x, as ~ x^''^^ - a manifestation 
of good quality of the Parisi- Toulouse scaling approximation. 



as a good approximation at finite x. The correction to 
scaling functions Sfh(z) and Sp{z) are shown in Fig. 

Numerical solution of the RSB equations. To complete 
the proof we have to demonstrate that the exact solution 
of Parisi equations indeed flows to the fixed point, rep- 
resented by our scaling solution, which we had shown 
to be attractive. A simple way to do this is to solve 
the exact RSB equations [Eqs. (|5I6|) ] numerically on the 
interval x £ [xq , 1] , using the scaling solution p as the 
initial conditions at the point — s- 0. If the function 
c{x) computed in that way is consistent with our scaling 
ansatz, that is c{x) is approaching, for small x, the con- 
stant c = 0.410802... (and fh is approaching its scaling 
counterpart), then one can claim that the scaling region 
does exist. A similar self-consistency check could be done 
on the interval x € [0,xq], though it is not as simple. In 
addition to proving our point, the outlined procedure will 
allow us to compute certain low energy quantities char- 
acterized by a finite time scale x — 1. 

In practice, when integrating Eqs. (|5I6|) numeri- 
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FIG. 3: Function c{x) obtained by solving the Parisi equations 
numerically at T = 0. It quickly approaches the constant 
value c = 0.410802... as i ^ 0. The inset shows In [c — c{x)] 
plotted vs In x (circles) along with the linear analysis predic- 
tion 5.41 Ins + const (solid line). Notice that the function 
c — c{x) spans more than three orders of magnitude. 

cally, it is convenient to switch to a new variable 
t = lux and avoid numerical differentiation by com- 
puting c{x) in successive iterations as Ci-^i{x) = 
Ci{x) J dzp{x, z){m' {x, z))'^ . An initial guess for c{x) can 
be a rather arbitrary positive function. The numerical 
solution quickly converges at small x, while converging 
logarithmically slow for x above and near the breaking 
point Xmax The process can be speeded up by choos- 
ing judicially |l4| the initial guess of c{x), taking in con- 
sideration that it vanishes for x > Xmax- The results are 
shown in Fig. (j^J. The numerical solution reproduces the 
exact value of c to five digits. In the inset we produce a 
log- log plot of c— c(x) together with our stability analysis 
prediction of const • x^''^^ . These results demonstrate the 
excellent approximating quality, as well as the validity, 
of the scaling ansatz and its correction. 

Using the numerical solution of the Parisi equations, 
namely the functions p(l, z) and m(l, z), we compute the 
leading low temperature behavior of the zero field cooled 
susceptibility xo — /3[1 — ?(1)] 

/3xo = j dzp{l,z)[l-m\l,z% (13) 

and entropy S: 

+ j dzp(l,z) [ln(2coshz) - ztanhz] . (14) 

We found /3xo ^ 1.592±0.002 and (3^3 = 0.7210±0.0002, 
in agreement with known estimates [l^ . and, to our 
best knowledge, with precision exceeding all other meth- 
ods. These values were computed for demonstrational 
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FIG. 4: Universal distribution functions p{z) (instantaneous 
fields, bottom graph) and p{l,z) (frozen fields, top graph) 
obtained by numerically solving the Parisi equations at T = 0. 
The instantaneous fields are enhanced, relative to the frozen 
fields, by the value of the Onsager term /3xo at large z. 



We would like to note that the idea of existence of a 
scaling solution of some sort for the SK model is not new, 
one of the recent attempts was reported in Ref 
our paper however, for the first time, a proposed scaling 
ansatz is proven to become asymptotically exact in a cer- 
tain limit. The presented method is not limited to the 
SK model only. It is expected to be applicable and very 
useful for a wide variety of spin glass models which admit 
full RSB scenario, such as the Ising p-spin model or a re- 
cently introduced mean field description of the Coulomb 
glass 0,0,0]. The scaling ideas may also be promising 
for understanding quantum glasses, which are thought to 
be relevant to the most challenging problems of strongly 
correlated systems [ioj . 

The author acknowledges extremely enlightening dis- 
cussions with Vladimir Dobrosavljevic and Markus 
Mueller. 



purposes only, and with some effort a far better accu- 
racy can be achieved. We also estimated the breaking 
point value to be Xmax — 0.56 ± 0.02 which is in agree- 
ment with the literature. Finally we compute a uni- 
versal distribution function of the instantaneous fields 
p{z) = P(zT)/T\t~+o, where P{y) is experimentally mea- 
surable distribution of instantaneous fields corre- 
sponding to the tunneling density of states in the context 
of charge glasses: 

f dz' 

' J ^2^0 ^ ' 

coshz r (z-z'f /3xol 

X — — 7 exp — — . (15) 

coshz' ^ [ 2/3x0 2 J ^ ' 

The distribution function p{z) is plotted in Fig. to- 
gether with p{\,z). At large z the distribution of instan- 
taneous fields is shifted, relative to the distribution of 
frozen fields due to effect of the Onsager reaction term 
[11 m by 7/3x0 « 0.479. 

Conclusion. To summarize, in this paper we presented 
a scaling ansatz for the zero temperature solution of the 
SK model. We assessed the stability of the scaling solu- 
tion, both analytically and numerically, and found that 
it becomes exact in the limit of zero temperature and 
low energy, for T <C x ^ and y <C 1. Our ansatz 
enabled us to compute the slope of the gap in the local 
fields distribution function numerically exactly, that is 
up to an arbitrary requested precision. Using the scaling 
solution as the initial conditions for the RSB equations 
at a; — !■ 0, we could compute, with high precision, the 
asymptotic behavior of finite time scale quantities, such 
as the entropy and zero field cooled susceptibility. 
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